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Abstract 

We work out two different analytical methods for calculating the boundary of the Mott-insulator- 
superfluid (MI-SF) quantum phase transition for scalar bosons in cubic optical lattices of arbitrary 
dimension at zero temperature which improve upon the seminal mean-field result. The first one is 
a variational method, which is inspired by variational perturbation theory, whereas the second one 
is based on the field-theoretic concept of effective potential. Within both analytical approaches we 
achieve a considerable improvement of the location of the MI-SF quantum phase transition for the 
first Mott lobe in excellent agreement with recent numerical results from Quantum Monte-Carlo 
simulations in two and three dimensions. Thus, our analytical results for the whole quantum phase 
diagram can be regarded as being essentially exact for all practical purposes. 
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I. INTRODUCTION 



During the last few years, experiments on trapped dilute ultracold quantum gases led to 
the observation of many new novel properties of these quantum systems . Among 

these systems are ultracold bosonic gases trapped in the periodic potential of optical lattices. 
These led to a whole plethora of experimental possibilities on many-particle quantum physics 
as they represent model systems for solid-state physics with a yet unprecedented level of 
control 5, (], 3, Is, [9], 1(J U, Q, Q]- Optical lattices can be formed by using electromagnetic 
standing waves orthogonally aligned to each other, with their crossing point positioned at 
the center of a Bose-Einstein condensate. In this way, they generate a periodic potential 
where the atoms can move from one lattice site to the next due to the quantum mechanical 
tunnelling effect [4]. Because of the periodicity of the optical lattice, the single-particle 
energy spectrum has a band structure. At low enough temperatures, the thermal fluctuations 
are too weak to excite the atoms beyond the lowest band. Thus, this system can be well 
described by the Bose-Hubbard model 
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Such bosonic gases in optical lattices can exist in two different phases which can be 
chosen by tuning the depth of the potential wells generated by the optical waves. The Bose- 
Hubbard model states that the existence of these two phases is determined by the balance 
between the atom-atom on-site interaction U and the hopping amplitude t. When the on-site 
interaction is small compared to the hopping amplitude, the ground state is superfTuid (SF), 
as the bosons are delocalized and phase coherent over the whole lattice. In the opposite 
limit, where the on-site interaction dominates over the hopping term, the ground state is 
a Mott insulator (MI), as each boson is trapped in one of the respective potential minima. 
These different phases are observable, for instance, in time-of-flight absorption pictures 
which are taken after switching off the lattice potential. While the superfluid phase yields 
distinct Bragg-like interference peaks, the Mott phase is characterized by a broad diffusive 
interference pattern 0, Q] . 

A common approximation for calculating the SF-MI phase boundary uses a mean-field 
theory where the non-local Bose-Hubbard Hamiltonian is substituted by an effective local 
one 14j . An alternative way to obtain the SF-MI phase boundary at T = is based on 



a strong-coupling expansion as is worked out in detail in Ref. 



181 ] . However, compar- 



ing both analytical methods with recent high-precision Monte-Carlo data 19], as shown in 



2 



t/u 

0.04 
.03 



0.2 0.4 0.6 0.8 



FIG. 1: (Color online) Quantum phase diagram of the first MI-SF lobe at T = for the three- 
dimensional case. Dot-dashed blue line is the mean-field result [l4], dotted black line is from the 
third-order strong-coupling expansion [la ), and red dots are recent high-precision Monte-Carlo data 

m 

Fig. Q] for the three-dimensional case, we observe that the mean-field theory underestimates 
the location of the quantum phase transition, while the strong-coupling approach overes- 
timates it. Thus, in view of a more quantitative comparison with experimental results, it 



becomes indispensable to 
stance, Refs. |2Q 



21. 
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urther develop analytical approximation methods (see, for in- 
, 2^1). In particular, obtaining accurate analytical results for the 
phase boundary at arbitrary dimension d and lobe number n would yield new insight be- 
yond the purely numerical data provided by Monte-Carlo simulation. Furthermore, it would 
be favorable to develop analytical approaches which allow, at least in principle, systematic 
improvements to higher orders. 

To this end we present here two alternative analytical methods to approximately solve 
the homogeneous Bose-Hubbard Hamiltonian and calculate the properties of the zero- 
temperature quantum phase diagram, in particular the location of the MI-SF phase bound- 
ary. Both our approaches are technically based on a systematic expansion with respect to the 
hopping parameter t. This perturbative procedure is justified as the phase boundary occurs 
in three dimensions for small values of t/U (see Fig. CD). It turns out that this i-expansion 
can be worked out analytically up to higher orders which contain non-trivial information 
about the dimension d, and, therefore, improve considerably the mean-field result for lower 
dimensions. Note that our approaches essentially differ from the promising, recently devel- 
oped Bosonic Dynamical Mean-Field Theory (BDMFT) which is based on a systematic l/d 
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expansion [21|. Although BDMFT has the virtue of being nonperturbative in the system 
parameters t and U, its self-consistency equations can only be solved numerically. 

In Section II we briefly review the usual mean-field theory which already gives a good 
qualitative description of the MI-SF quantum phase transition. In Section III we work out 
our first method where we perform systematic variational corrections to the previous mean- 
field calculations by considering a different interpretation of the order parameter ip. Unlike in 
mean-field theory, where ip is determined from self-consistency relations, here tp is regarded as 



a variational parameter like in variational perturbation theory |25l . l26j . In this new approach, 
the order parameter ip is found by extremizing the grand-canonical free energy, although 
self-consistency relations no longer apply in general. Afterwards, we present our second 
approach in Section IV which is based on a standard field-theoretical method. Here we 
add spatially and temporally global source terms to the usual Bose-Hubbard Hamiltonian 
in order to break the global U(l) symmetry. The Legendre transformation of the grand- 
canonical free energy obtained from this modified Hamiltonian defines an effective potential 
which is used for determining the MI-SF phase-boundary. Both analytical approaches are 
calculated up to second order in the hopping parameter t in Section V, while detailed 
technical calculations are relegated to the Appendices. Finally, we present in Section VI the 
resulting improved quantum phase diagrams and compare them with previous findings. In 
particular, we show that our analytical results have an accuracy for the first MI-SF lobe in 
two and three dimensions which is comparable with Monte-Carlo data. Therefore, we are 
confident that our analytical results are essentially exact for any Mott lobe in more than 
one dimension. 

II. MEAN-FIELD THEORY 



In the present work we deal with a system of spinless bosons in a homogeneous infinite 
cubic lattice of arbitrary dimension d, which can be well described by the Bose-Hubbard 
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Hamiltonian [JJ, [JJ, [16|, [17 1 



H B}1 = -tJ^aldj + Ho (1) 
with running over the nearest neighbor sites and with the on-site Hamiltonian 



Here hi = aja, represents the number operator at site i, and /i denotes the chemical potential. 
Furthermore, t is the hopping energy which characterizes the tunneling of an atom from one 
lattice site to a neighboring one. Furthermore, U denotes the on-site energy which describes 
the strength of the interaction between two atoms at a given lattice site. When the depth 
of the lattice wells is increased, the hopping energy t decays exponentially fast, whereas the 



on-site energy U increases algebraically 



27]. 



A standard approach to approximately solve the Hamiltonian ([I]) uses a mean-field ansatz 
[yj]. To this end, the non-local hopping term in (CQ) is substituted by a sum of single-site 
terms, thus yielding the following mean-field Hamiltonian: 

H MF = -t(2d) (V ^ + #i - ^) + H . (3) 

i 

The thermodynamical quantities are then calculated from the grand-canonical free energy 

F MF (r^) = ~^Z MF (r,^P), (4) 
where the grand-canonical partition function is given by 

Z MF (f, V) = Tr [e-^^'^j . (5) 

The additional complex parameters ip and ip*, which effectively describe the influence of the 
neighboring sites, have to be self-consistently determined according to ip* = (a\), ip = (Sj). 
Note that imposing these self-consistency conditions is equivalent to extremizing the grand- 
canonical free energy with respect to ip and ip* : 

dF MF 



dip 
dF MF 









(6) 



dip* 

Near the phase boundary the order parameter ip is small, so the grand-canonical free energy 
can be Taylor expanded with respect to the order parameter in form of a Landau expansion 



2a M: 



FMP,1>) = N s K F (T) + af F (T)\iP\ 2 + af F (T)\iP\* + •••]. (7) 

Here N s is the total number of lattice sites within the system. If 04 (T) > and 02 (T) 
changes its sign depending on the values of the respective parameters t, U, and /1, then the 
system exhibits a second-order phase transition. The phase boundary is given by points in 



the parameter s pac e where a^ F (T) = is valid. At zero temperature the resulting phase 
boundary reads [141 ]: 

Note that this result of mean- field theory becomes exact in the limit d — > oo 

III. VARIATIONAL METHOD 

As the mean-field results differ considerably from the latest Quantum Monte-Carlo results 
(see Fig. [I]), it is necessary to develop new analytical approaches for studying bosons in 
optical lattices. Therefore, the objective of the present work is to generalize the mean-field 
approach and to calculate systematic corrections which are due to quantum fluctuations. 
Our first method is based on introducing an artificial smallness parameter rj in the Bose- 
Hubbard Hamiltonian according to 

H (77) = H MF + V (h bu - H MF ) , (9) 

which can be explicitly written as: 

H{r], r, ^) = H -tr)J2 ^ - 2dt ^ ty a * + - M 2 ) • (10) 

Note the limiting cases rj = and 77 = 1, where the Hamiltonian fflQj) reduces to H (rj — 0) = 
Hmf and H {rj — 1) = -^bh, respectively. Furthermore, the order parameter ip is treated as 



a variational parameter like in variational perturbation theory (VPT) [25|, |26j. Using we 
perform a Taylor expansion up to the iVth order of the grand-canonical free energy in r) and 
obtain, as in the mean-field case, a Landau expansion for this iVth order grand-canonical 
free energy: 



F^( V ,r^) = N s ai N \7 1 ) + ai N \r ] M 2 + a{ N \^ + ... . (11) 



Here is the truncated expansion in r/ up to order N of a 2p (ri) which is given by the 
expressions 

\ 2 ( V ) = {2dt)\l - V f E^=o(-^) m «2 m) + 2A(1 - v), 

a 2p ( v ) = ^dt) 2p a-v) 2p E^ =0 (-tv) m cy { 2 ; ) ; p^l. 



When we set p = 1 and rj = 1, the truncated expansions reduce to 
aP(r] = l) = af\2d)H 2 + 2dt, 



a i 2 N) (r ] = l) = (-l) N e 



4 N k N + of" V- 1 



(13) 



iV > 1. 



Finally, we find the phase boundary at the points where we have a^\ri = 1) = 0. Thus, we 
conclude from f[T3"j) that the phase boundary is given by 



I c 



(2rf)4 0) ' (u) 

,, n (iv-i) v ; 

f( N ) - _^2_ _ . AT > 1 

r c — (lv)- i iV ^ 1. 
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Before we calculate explicitly the respective perturbative coefficients c4 m (ED; we intro- 
duce in the next section our second method for determining the location of the quantum 
phase transition. 

IV. FIELD-THEORETIC METHOD 

The second method developed here in order to improve the analytical results for bosons in 
optical lattices is not based on any mean-field theory. Instead of this, as a starting point, we 
consider the Bose-Hubbard Hamiltonian with additional source terms, which are spatially 
and temporally global: 

h br (j\ j) = -tJ2 a l a * + ( J * a « + Ja + (is) 

(m) * 

The grand-canonical free energy is then calculated in a power series of both the hopping 
parameter t and the sources J, J*. This leads to 



F(J*,J,t) =N S ^F (t) + J2c 2p (t)\J\ 2p j (16) 
with the expansion coefficients 

oo 

C2 P (t) = ^(-t)"«g } . (17) 



n=0 



We observe that due to the similarities between the Hamiltonians (fTUl) and (|T5l) . the coeffi- 
cients a^p which appear in the expression ffTTj) are identical to the ones in Eqs. (1121) . 
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In this approach we define the order parameter ip, as usual in field theory 25|, l28j . 
according to 

. . i dF{r,j) t 1 dF(j*,j) 

A subsequent Legendre transform of the grand-canonical free energy yields the effective 
potential 

r(r,^) = F/N B -if)*J-^J*- (19) 

Thus, the external sources can be written as derivatives of the effective potential 

^ J , J*. (20) 



dip* ' dip 

From Eq. (f20l) we read off that the physical limit of vanishing currents is obtained by simply 
extremizing the effective potential with respect to ip and ip*: 

W = • l = ' (21) 

Using (IT61) . the effective potential (IT9|) is written as a power series of \ip\ 2 : 

r(r, ^ t) = F Q (t) - * + + ■ ■ ■ . ( 22 ) 

C 2W c 2 ( i tj 

The respective coefficients can be determined as a power series in t. For instance, the 
coefficient of \ip\ 2 turns out to have the following hopping expansion: 

„(2)' 



1 + -L-t + 



c 2 (t) J® 




(o) 
a 2 



r + • • • > . (23) 



l 2 I "2 

This expansion of l/c 2 in power series of t is equivalent to a resummation of c 2 such that 
it has a divergency at the phase boundary. The presence of such a divergency comes from 
the long-range correlations of the system and is an essential feature for the occurrence of 
the phase transition. Thus, the quantum phase boundary is found by setting l/c2(t c ) = 0. 
This procedure gives us an algebraic equation in t c whose degree depends on the order of the 
expansion in t. Such an algebraic equation can have, in principle, many real roots depending 
on its degree. However, only the smallest root must be considered as physical. The other 
roots must be discarded as they are artificially introduced within our present method whose 
validity is restricted to small values of t. Following these criteria we find in first hopping 
order: 

(o) 

(?» = (24) 



a 



2 
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For the second order we get correspondingly: 

( ' 2) = 2(i^f) + W^W) ^ 52) ' (25) 

where we have introduced the reduced quantities a\ = a^/a^ and «2 — ot f/af. Note 
that only the smallest root was taken into account in (1251) . 

V. PERTURB ATIVE CALCULATIONS 

In this paper we restrict ourselves to work out the zero-temperature limit of our two 
analytical approaches. Therefore, we can use the standard Rayleigh-Schrodinger perturba- 
tion theory in order to obtain the ground-state energy corresponding to the Hamiltonians 
presented here. Obviously, the Hamiltonian ffl5l) is converted into the Hamiltonian fflUj) by 
the following transformation of its parameters: 

t ► 7]t, 

J — ► 2dt(r] - l)ip, 

J* — ► 2dt(r] - 1)V*, (26) 

H — ► H + 2dt(l -r])\i)\ 2 . 

This fact enables us to perform all calculations using the Hamiltonian of the field-theoretic 
method ffl5l) and then obtain the results corresponding to the variational method by re- 
defining its parameters. Therefore, all we need to calculate is the ground-state energy of the 
Hamiltonian ffl5l) . 

As the ground-state energy of ( 1151) is an extensive quantity, its calculation as a power 
series in the hopping parameter t can be performed by applying the linked cluster method 
(see Appendix As can be seen in f lTTl) . the coefficients are the fundamental blocks for 
calculating all quantities which are related to the ground-state energy. In our version of the 
linked-cluster expansion each of the coefficients ot$ can be represented by a set of diagrams. 
Such diagrams are composed of oriented lines linked at their ends to points which represent 
the vertices of the lattice. The diagrams are embedded in the lattice and, therefore, their 
topologies are defined by the topology of the underlying lattice. The linked-cluster theorem 
states that only connected diagrams contribute to extensive quantities like the ground-state 
energy. 
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As explained in detail in Appendix [Bl for a given diagram, there are lines which have 
both ends linked to neighbor points of the lattice (internal lines) and there are lines which 
are linked to a point of the lattice by only one of its ends (external lines). The internal lines 
are associated with the hopping of the atoms between two neighboring points and, therefore, 
the order in the hopping parameter t of a diagram is given by the number of its internal 
lines. The external lines can be directed into or out of a lattice point. If the line is directed 
into the point it represents a creation operator acting at this point which is the coefficient 
of J in ffT5l) . In the opposite case where the line is directed out of the point it represents an 
annihilation operator acting at this point which is the coefficient of J* in ( |T5l) . Analogously 
to the internal lines the power of J is given by the number of lines going into the diagram 
while the power of J* is given by the number of lines coming out of the diagrams. As the 
ground-state energy of ffT5j) depends only on J* J, the number of ingoing lines must be equal 
to the number of outgoing lines. 

(n) 

The location of the MI-SF phase boundary is determined by the coefficients a 2 whose di- 
agrams involve two external lines. The zeroth-order coefficient corresponds to diagrams 
with no internal lines, thus, for topological reasons, only one diagram can contribute: 

af= » • » • (27) 

In first order the coefficient corresponds to one internal line, so also there only one 
diagram is possible: 

a W = (2d) (28) 

Here 2d is the corresponding lattice number which can be explained as follows. Once the 
first vertex is chosen somewhere in the lattice, there are precisely 2d possibilities for the 
second vertex. The situation is more complicated for the second-order coefficient , as in 
that case two topologically different diagrams contribute: 



oj 2) = (2d) (2d- 1) » »+(2rf) » V » (29) 

In Appendix Owe develop a second kind of diagrammatics which is useful for simplifying 
the calculations of the coefficients (I27]) - (I29~]) within time-independent perturbation theory. 
This explicit evaluation is relegated to Appendix [D] which yields the following expressions: 
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FIG. 2: (Color online) Quantum phase diagram of the first MI-SF lobe (n = 1) at T = 0. Solid 



green lines are results from our field-theoretic method, dashed purple lines are results from our 
variational method, dot-dashed blue lines are from mean-field theory dotted black lines are 



from the third-order strong-coupling expansion 



la ], and red triangles are the numerical data. For 



the one-dimensional case the numerical data stem from Density-Matrix-Renormalization Group 
e the data for two and three dimensions are obtained from Quantum Monte- 



calculations 



29|], whi 



Carlo simulations 



19 



a 



(o) 



U(b-n)(b+l -n)' 



(30a) 



a 



(i) 



2(2(6 + 1) 



U 2 (b-n) 2 (b+ 1 -nf 



(30b) 



a 



(2) 



2d {2d(b + lf{b - 2 - n)(b + 3 - n) + n(b - n)(b + 1 - n)(l + n) 



x (4 + 36 + 2n) [-3 - 2n + 2{b 2 + b - 2bn + n 2 )] } 

/ [U 3 (b -n-2)(b- nf{b + 1 - nf{b + 3 - n)] , (30c) 
where we have used again the abbreviation b = fi/U. 



VI. QUANTUM PHASE DIAGRAMS 

In order to calculate the quantum phase diagram by using the variational method, we 
have to substitute the coefficients fl30l into Eqs. f|T^|) . while for the field-theoretic method 
they are substituted into Eqs. (1241) and (|25|) . 
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TABLE I: Relative deviation of our analytical findings from the Quantum Monte-Carlo data in the 
position of the first MI-SF lobe tip in two and three dimensions. 





mean-field theory 


strong coupling 


variational method 


field-theoretic method 


d = 2 


28% 


13% 


13% 


6% 


d = 3 


16% 


24% 


5% 


3% 



In the variational method, the zeroth-order phase boundary reads 



tf) ~ (2d)«r (31) 

whereas the field-theoretic method yields in lowest order 

(0) 

t? = (32) 

By applying the identities (1B2I) we see that Eqs. (l3Tj) and (1321) are actually one and the same 
expression, which reduces with the help of (130 aft to the mean-field result (jSJ). 

The first non-trivial quantum correction to the quantum phase boundary in both ap- 
proaches becomes available just at the second-order level, where the formulas (|14|) and (125]) 
have to be applied together with (1301) . In Fig. [2] we compare, for the first MI-SF lobe, our 
two methods with the mean- field phase boundary , with the third-order strong-coupling 



expansion [18( , and with recent numerical data obtained by using the Density-Matrix Renor- 



malization Group technique 31] for one dimension as well as the Quantum Monte-Carlo 



30] simulations for two and three dimensions. In Tabled] we compare our two analytical 
methods quantitatively with the Quantum Monte-Carlo data for the position of the first 
MI-SF lobe tip for two and three dimensions. Thus, we can conclude from Fig. [2] and Table 
HI that our analytical methods have provided very accurate results for the first MI-SF lobe 
in comparison with Monte-Carlo data. Therefore, we expect that they are essentially exact 
for any MI-SF lobe in more than one dimension where no Monte-Carlo simulations have 
been yet performed systematically. In addition, we read off from Fig. [2] and Table [I] that 
our field-theoretic method turns out to give better results for the MI-SF phase boundary 
in two and three dimensions when compared with the variational method. Despite of this 
discrepancy between our two analytical approaches, their theoretical formulations suggest 
that both converge to the true result once higher hopping orders are taken into account. 
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FIG. 3: (Color online) Phase diagram of MI-SF lobes at T = for different dimensions. Solid green 
line are results from our field-theoretic method, dashed purple lines are results from our variational 
method, dot-dashed blue lines are from mean-field theory, and dotted black lines are from the third- 

la |. The non-physical 10-dimensional case is included to show 



order strong coupling-expansion 
that our two methods converge to the mean-field theory 
general grounds [lfj]. 



14l ] in the limit d — ► oo, as is expected on 



For d — 1 the quantum phase boundary of the Bose-Hubbard model is more complicated 



as it is a Kosterlitz-Thouless type of phase transition If], 32]. This non-analytic behavior is 
reproduced in Fig. [2] quite well by both the precise Density-Matrix- Renormalization Group 
results and the strong-coupling expansion (see also Ref. 33|). However, our two analytical 
methods cannot deal with this. Our variational method yields, at least, a continuous phase 
boundary, but our field-theoretical method leads to a finite interval of the chemical potential 
where no real solution for the phase boundary exists. This finding is insofar consistent as 
both methods are expected to be applicable only for small values of the hopping parameter 
t. 
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In Fig. [3] we can clearly observe how the phase boundaries of our two methods approach 
the mean-field phase boundary as the dimension d of the system increases which is consistent 
with the fact that at large dimensions the mean- field theory becomes exact 16|]. Such an 
agreement with the mean-field theory for d —>■ oo can also be checked directly in Eqs. (|14l) . 
(|24p . and (1251) by using the explicit results for 2 , and c£p from Eqs. (l3"Uj) . 

Finally, we observe that so far our two methods yield a phase boundary for d = 3 
dimensions which turns out to be analytical at the lobe tip. This finding is consistent 
with the theory of critical phenomena as a quantum phase transition in d spatial dimensions 



belongs effectively to the universality class of a standard phase transition in d+1 dimensions 
I4I \m . In contrast to that, the strong-coupling expansion of Ref. 18] leads in each 
order to a pronounced artificial cusp at the lobe tip. At present, it remains open to resolve 



the analytical structure of lobe tips via QMC simulations [19j. This is certainly demanding 
as the lobe tip is the most sensible region of the Mott lobe with respect to finite-size scalings. 
We expect that also at higher orders our two methods will show an analytical lobe tip. This 
would make them then ideal tools for dealing with experimental situations where the MI-SF 
phase boundary is crossed at a fixed particle number per site. 



VII. CONCLUSIONS AND OUTLOOK 



In this work we have presented two alternative analytical methods which improve sub- 
stantially the quantitative results for the MI-SF quantum phase boundary for dimensions 
larger than one as shown in Fig. [2] and Table [fl For the first Mott lobe our results are in the 
immediate vicinity of recent high-precision Monte-Carlo simulations. Therefore, it would be 
interesting to investigate with future Monte-Carlo simulations whether the accuracy of our 
results increases or decreases for higher Mott lobes. 

We expect that both methods converge to one and the same result if higher orders in the 
hopping energy would be taken into account. However, as can be seen in Table HI at least 
up to the second order in the hooping parameter t, the field-theoretical method turns out 
to converge faster than the variational method. Another advantage of the second method 
is that the order parameter ip is, by definition, the square root of the condensate density, 
while for the first method the physical interpretation of ip remains unclear. 

In addition, we remark that both theories in Sections III and IV have been formulated 
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so general that they are also applicable for finite temperatures. Thus, it should soon be 
possible to improve the finite-temperature mean-field MI-SF phase boundary which was 



recently derived in Refs. 



34 



35 



36|. 



Both approaches presented here for a homogeneous optical lattice can be extended in a 
straight-forward way to the experimental situation where an additional harmonic confine- 
ment potential is superimposed to the periodic potential. Within a Thomas-Fermi approxi- 
mation the overall harmonic potential is taken into account by introducing a local chemical 
potential. Thus, the particle density distribution of such an inhomogeneous system follows 
from cutting the Mott lobes of a homogeneous quantum phase diagram with a horizontal 
line. This yields to a wedding cake structure which consists of alternating concentric insu- 
lating and superfluid layers . A precise description of the respective layer widths is 
of fundamental importance in view of a quantitative analysis of time-of-flight data j^ . 

Finally, we conclude with the observation that our present approaches are so far restricted 
to the Mott-insulating regime. In order to extend them to the superfluid regime, where we 
have a non-vanishing order parameter, would necessitate to determine higher order terms in 
the Landau expansion. In addition, to have access also to local statistical quantities as, for 
instance, correlations functions it is indispensable to develop a Ginzburg-Landau expansion 
where the spatially and temporally homogeneous order parameter is generalized to an order 
parameter field [381]. 
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APPENDIX A: LINKED-CLUSTER THEOREM 



Different versions of the linked cluster expansion have so far been applied to both classical 
and quantum lattices systems. This section is dedicated to applying the linked-cluster 
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theorem [39] to the Bose-Hubbard Hamiltonian with currents. Our first step is to introduce 
the more general Hamiltonian: 

Hbh = -y^%,j>alaj + ^2 ("o^Wi^i ~ ^laA + ^ (j?ai + Jihfj . (Al) 

The difference between the above Hamiltonian and the Hamiltonian (j!5p resides in the more 
general coefficients J*, and Jj. Now the hopping parameter t can have different values 
for each pair of nearest neighbor lattice sites and the sources can have different values 
for each lattice site. 

Now imagine a given extensive quantity G associated with the system described by (1A11) . 
Such a quantity can be expanded in a power series of the indices 

{"(i, 3 )} (*»i> 

Here {^(ij)} denotes the set of all coefficients associated with the nearest neighbor pairs and 
{ n (i,j)} is a se t of discrete indices also associated with the nearest neighbor pairs with each 
n/ij\ running over the nonnegative integers. Without loss of generality, we can assume that 
Q = if tfaj) = 0, which is equivalent to assume g{o,o,-} — 0. 

A cluster is defined as any nonempty set of neighbor pairs . With this definition f]A2|) 
can be written as 

Q({t {t , j) }) = J2w(C), (A3) 

c 

where the sum runs over all possible clusters C. The cluster weight W(C) contains all terms 
in (IA2I) which have at least one power of ta j\ for all contained in C and no powers 
of any other tujy The linked cluster theorem states that the weight W(C^) associated with 
a disconnected cluster d is zero and, therefore, only connected clusters must be taken 
into account in f]A3|) . Here a disconnected cluster represents any cluster which is just a 
union of disjoint nonempty subclusters. In order to verify the linked cluster theorem we 
can suppose that only the coefficients t^) associated with the pairs contained in two 
disjoint nonempty clusters C\ and C 2 are nonzero and all other t^j) are zero. The additivity 
of Q then implies that 

Q(C 1 UC 2 )=Q(C 1 )+Q(C 2 ). (A4) 

It means that in (1A2|) there is no term which includes t^j) for all contained in C\ UC 2 , 
therefore we have W{C\ U C 2 ) = 0. 
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Although in (lAlj) a different t^j) is used for each pair in our applications we are 

interested in the special case where tuj) = t for all In this case clusters related 

to each other by symmetries such as translations, reflections, and rotations give the same 
contribution for Q. Each set of equivalent clusters can be represented by a graph G and 
the lattice constant L(G) is the number of topologically equivalent clusters represented by 
G per lattice site. Using these definitions Eq. (1A3[) reduces to 

Q(t) = N B J2L(G)W(G), (A5) 

G 

where N s denotes the number of lattice sites. 



APPENDIX B: DIAGRAMMATIC CLUSTER REPRESENTATION 

According to Appendix |X] a given cluster C is associated with a weight W(C) which is 
given by the terms in Eq. (1A2j) which have the product of coefficients tuj\ corresponding 
to all pairs in C, and without coefficients corresponding to pairs not contained in C. Each 
of these terms contributing to a W(C) can be further expanded in a power series of the 
currents Ji and J*. If we set now constant values to Ji, and J* and consider G as being 
the correction to the grand-canonical free energy due to t, J, and J*, we obtain the power 
series (TT6|) and (fPTl) . 

As we can see in flAlj) each coefficient tu j\ is associated with an annihilation operator dj 
and a creation operator a\, while each current Jj and J* is associated with a creation and 
annihilation operator a\ and di, respectively. Therefore, a set of diagrams composed of n 
internal lines and 2p external lines linked to points, which correspond to one and the same 
cluster, can be associated to each coefficient . With this each coefficient o^p can be 
written as a sum of diagrams where each diagram is multiplied by its lattice number. With 
this we obtain, for example, for n = and n = 1 the diagrams ([27]) and (|28|) . 

An interesting property of the diagram (|28|) is that it is one-particle reducible, and there- 
fore, similarly to the usual Feynman diagrams used in quantum field theory, it factorizes 



into its one-particle irreducible contributions 25|, |28|. Thus, the diagrams ([27]) and (|28|) are 
related via 



(Bl) 
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which yields in total 

<#> = (2d)(a^f. (B2) 

APPENDIX C: DIAGRAMMATIC NOTATION FOR TIME-INDEPENDENT 
PERTURBATION THEORY 

In this section a diagrammatic version of the usual Rayleigh-Schrodinger perturbation 
theory is presented. Using this approach the diagrams mentioned in the previous section 
are further decomposed into simpler contributions which are then easily computed. 



In the conventional time-independent perturbation theory 40| the aim is to solve the 
eigenvalue problem 

= (ci) 

where H = Hq + XV, in power series of the smallness parameter A: 

oo 

i*»> = E Ai i*®>> (° 2 ) 

i=0 
oo 

E n = J2 X ' E n- (C3) 

i=0 

Solving the Schrodinger equation (ICip with (1C2P and (103|) leads to the following recursion 
formula for % > 0: 



n \ n \" \ n 1 1 

l*n / 2^^ m > p (0) _ p(0) " ^ m ' p (0) _ p (0) • 

m^n j=l m^n ^« ^"i 

The first three terms of the eigenvalue expansion read explicitly: 
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EV> = (*W\V\^), (C5) 
= E ^^(o) ^ ^!^)^!^!^), (C6) 



n p (0) _ p (0) 771(D) p (0) \ ™ 1 1 ^m 2 /V^m 2 l " \*mi/\*mx\ v \*n I 

- E , F( o) ~ f( o), 2 (^ 0) i^i< 0) )(^ 0) i^i^!)(^!i^i^ 0) )- ( C7 ) 

mi^n (A« E mi ) 

The results in Eqs. ( 1C5l) -( JC7l) can be graphically represented by introducing a diagrammatic 
notation: 



= , (C8) 

E® = ^+ (C9) 



#i 3) = - (CIO) 



Here the diagrams in Eqs. flC8h — flCIO)) must be read from right to left and are interpreted 
according to the following rules: 

• Each dot represents an interaction V. Therefore, the number of dots contained in a 
diagram defines the respective power of A this diagram corresponds to. 

• The internal lines, which appear between two consecutive dots, are associated with 
the factor ^2 mj L n ^ (0) ^oyy \^m)(^m\, where q denotes the number of lines linking 
two given consecutive dots. 

• If we have more than one line between two consecutive dots, each extra line is associ- 
ated with an extra disconnected part of the diagram. Note that each diagram has a 
prefactor (— where s is the total number of its disconnected parts. 

• The external lines are associated with the unperturbed bra and ket (vP^I • • • |^ , n°' ) )- 



For example, the diagram (108j) has just one dot and the external lines, the diagram (JC9J) 
consists of two dots, one internal line, and the external lines. Equation ( 1C10I) contains two 
diagrams, where the first one has three dots with two internal lines and the external lines. 
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The second one represents a disconnected diagram where the lower part has two consecutive 
dots with two internal lines between them. The extra internal line is associated with the 
upper disconnected part which has just one dot without internal lines. Using the prescription 
above we can easily reconstruct the formulas (jC5p — (jC7j) . 

Following this prescription the i-th term En of the perturbative expansion of the energy 
E n can be constructed by summing all the diagrams with % dots and multiplying each diagram 
by a factor (— where s is the number of disconnected parts of the diagram. 

Thus, the 4th order term of the energy reads diagrammatically 



E^=-m-9-++ — _ t0t _ # _"_ # _ t0t _"_ i0t _ + • • ( cn ) 

Using the prescription to write down the explicit formula for E^ 1 we have: 

K (4) = V V V /»T' (0) I V 1 mi I- { mi 1 V 1 1112 ' { m2 1 V 1 m3 '1 ma 1 V\¥ 0) ) 

Z_> \ V p (0) _ p (0) p (0) _ p (0) p (0) _ p (0) l*r» / 



|vl> (0) \/vl/ (0) l |v|/( h/v|>( )| 

E E (*l 0) M^ 0) )(^ 0) r 1 ^ 7 (C12) 



|vl>(0)w^(0)| lvl>( )\/vl/ (0) l 

E E t 1 2 mi 0) )(*i 0) r 2 a)( y 

I "till -^w ) 

+ E (< 0) l^l^ 0) )(^°Wl< 0) )(^ 0) |V- T^^f^ ^l^)- 
roi^n ( £ mi — -En 



<< »<- << I - t;; j ^ ^ n 

r (0)w^(o) 

-(0) _ p (0) 
'mi "n 

Now we can construct diagrammatically any term of the perturbative series without using 
the recursion formulas (1C4I) explicitly. Note that we read off from (1040 a recursion formula 
for the total multiplicity M k of all diagrams in a given perturbative order: 

k-i 

M k+1 = M k + E M k _ l+l M h (C13) 

Starting from the initial value M\ = 1 corresponding to ( 1C8I) . its recursive solution gives 
M 2 = 1, M 3 = 2, and M 4 = 5, in agreement with flC9[) — (1C11H and yields furthermore 
M 5 = 14 and M 6 = 42 for the next two orders. 
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This approach can also be used if more interactions are present. For example, if we have 
H — Hq + XV + <jW, then we can introduce the two vertices 



V -> ♦ , (C14) 



W -> ♦ . (C15) 



With this we obtain for the perturbative expansion 

m 1 2 / 1 2 

E n = + A — • — + a— •— + Xa (■ 





2 1 

2 N 

+ •••. (C16) 

12 1 1 7 

APPENDIX D: EXPLICIT EVALUATION OF DIAGRAMS 

The purpose of this section is to evaluate the 'arrow' diagrams presented in Appendix |B] 
by using the 'line' diagrams presented in Appendix O 

As we saw in Appendix |B] each line in an arrow diagram represents a process which is 
associated with the operators Sj, a\, or aja^ with i and j being nearest neighbor lattice sites. 
In the line diagrams each of these processes are represented by a point. Let us take as an 
example the following arrow diagram: 

^3^2 ^ 

i j 

Here we labeled the processes with the numbers 1, 3, and 2 which, according to (ID 11) . 
are associated with the operators dj, ajdj, and dj, respectively. This diagram belongs to 
a cluster which consists only of two neighbouring sites i and j. Therefore, we define an 
effective Hamiltonian which takes into account only the relevant sites % and j, as well as the 
relevant interaction terms corresponding to oj, djdj, and aj as follows: 

H ea = Hi + ^ + V X + V 2 + V 3 . (D2) 
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Here the Hamiltonians Hi and Hj are defined according to (|2J) and have the same set of 
eigenvalues given by 

e n = ^(n 2 -n)-fin (D3) 

and the interactions read 

Vi = Ja\ , V 2 = J*aj , % = -ta\a 3 . (D4) 

The value of fIDll) is given by that term in the expansion of the ground-state energy corre- 
sponding to ( 1D2I) which is proportional to —tJ*J. According to our discussion in Appendix 
ICl this term can be split using line diagrams. Therefore the diagram (IDlj) can be written 
as follows: 

132 123 132 

= — 1 h more 4 permutations. (D5) 

Each line diagram can now be evaluated using the rules presented in Appendix [B] As an 
example we consider the line diagram 

13 2 

01,3,2 = • • • , (D6) 

which must be evaluated as follows: 

• Reading the diagram from right to left we find first the dot number two which is 
associated with the operator a,j. It annihilates one particle at site j generating the 
factors y/n and l/(e n — e n _i) and leaving n — 1 at this site: 

01,3,2 = Vn — — • (D7) 

• The next dot is number three which is associated with ajSj. Thus a particle is created 
at site j while another one is annihilated at site i. This leads to two factors \fn and 
the factor l/(e n — e n _i): 

r- 1 1 / x 

0i,3,2 = Vn — n — . (D8) 

• The last dot is number one which is associated with at. It creates the last particle at 
site % and yields the additional factor \/n. Thus, we obtain at the end: 



1 



01,3,2 = n 2 - -. (D9) 
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Now we present a complete list of all arrow diagrams up to the second hopping order 
which have been evaluated using line diagrams. In order to simplify the notation we will 
use the definitions 



A +p — e n — £n+p, A p — e n — e n _ p , (D10) 
A+ = A +1 , A~ = A -1 . (Dll) 



We start with the calculation of the arrow diagram (1271) which decomposes into two line 
diagrams: 

12 12 2 1 

fr • » = • — • + • — • (D12) 



Both contributions in (ID12I) lead to the following expressions 



thus yielding in total the result 



1 2 1 

• • = + (D13) 

2 11 , 

• • = n— , (D14) 

A 



n n + 1 .„ 

F + — < D15) 



Combining (IDlOj) and (1D11I) with (1D3j) and (1D15|) thus results in (I30ap with the abbreviation 
b = n/U. 

Correspondingly, the arrow diagram in (|28|) factorizes according to (1B1I) . so we obtain 
with (TDT51) 

n n + 1 



,- + — ) • < D16 > 

which leads to the result fl30bl) . 

The first arrow diagram in (|29j) factorizes in a similar way, yielding 

i n + 1 \ 3 . 

x + — ) • < D17 > 

Now we return to the second arrow diagram in (1291) which decomposes in line diagrams 
according to: 



4 A3 • • 

V » >♦> ! 



- (D18) 
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Here the non-numbered dot diagrams represent a sum over all numbered diagrams with the 
respective topology. Their explicit evaluation yields 



1 3 2 4 n(n + l) 2 

# * * * (A+) 2 (A+ + A-)' 
1 4 2 3 n(n + l)(n + 2) 



* " A+(A+ 2 + A-)(A+ + A-)' 
1 3 4 2 (n + 1) 3 



(A+) 



3 



3 14 2 n(n + l) 2 



' (A+) 2 (A+ + A~)' 
4 13 2 n(n + l)(n + 2) 



A+(A+ 2 + A-)(A+ + A-)' 



2 4 13 n 2 (n+ 1) 



(A-) 2 (A+ + A-)' 



(D19) 
(D20) 
(D21) 



HI 2 . - n(n + l)(n + 2) 

" (A+) 2 (A+ 2 + A-) ' [UZZ) 
3 12 4 n 2 (n+ 1) 

™- = A +( A ++ A-r (D23) 

till - n(n + l)(n + 2) 

_ (A+ 2 + A-)(A+ + A-) 2 ' [UM) 



(D25) 
(D26) 



2 g j a - "(n-l)(n+l) (D27) 
_ A-(A+ + A- 2 )(A+ + A~)' [UZ{) 



(D28) 



32 14 n(n-l)(n + l) 

" (A+ + A- 2 )(A+ + A-)^ 

4 2 13 l) 2 
a a a a — = — 

A-(A+ + A-) 2 ' 

2 3 4 1 n (n- l)(n+ 1) 
# # # a — = — — - 

(A-) 2 (A+ + A- 2 )' 
2 4 3 1 n 3 

= (A 3 ) 3 ' 



" ii = „ + : w2 ,;u , ;,, 2 , (D29) 

(D30) 
(D31) 
(D32) 



3 2 4 1 n(n-l)(n + l) , 

= F^WTF) - (D33) 



4 2 3 1 n 2 (n + 1) 



(A-) 2 (A+ + A-)' 



(D34) 



4 2 3 1 n 2 (n + l) _ or , 

= (A-) 2 (A + + A-) ' (D35) 
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3 


4 


n 2 (n + l) 
(A ) 2 (A+ + A ) 




•- 

1 


-• — 

2 


(D36) 


4 


3 


n 2 (n + l) 
(A ) 2 (A+ + A ) 




• 


• — 

2 


(D37) 


3 


4 


n(n + l) 2 
(A+)^(A+ + A ) 




— • 

2 


• — 

1 


(D38) 


4 


3 


n(n + l) 2 
(A+) 2 (A+ + A-) 




•- 

2 


-• — 

1 


(D39) 


1 


2 


n 2 (n + l) 
A~(A+ + A~) 2 




•- 

3 


-• — 

4 


(D40) 


2 


1 


n(n + l) 2 

\ I / \ I i \ \o ! 

A+(A+ + A~) 2 




• 

3 


• — 

4 


(D41) 


1 


2 


n 2 (n + 1) 
. A~(A+ + A-) 2 ' 




•- 

A 


-• — 

Q 
O 


(D42) 


2 


1 


n(n + l) 2 
A+(A+ + A-) 2 ' 




•- 

4 


-• 

3 


(D43) 



Thus, combining flDTBD with flDTOl) . flDTTD and flDT8l) - flD43|) finally yields (I50c|) . 
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